Experimental measurement and numerical modeling of deformation behavior of breast cancer cells passing through constricted microfluidic channels

During the multistep process of metastasis, cancer cells encounter various mechanical forces which make them deform drastically. Developing accurate in-silico models, capable of simulating the interactions between the mechanical forces and highly deformable cancer cells, can pave the way for the development of novel diagnostic and predictive methods for metastatic progression. Spring-network models of cancer cell, empowered by our recently proposed identification approach, promises a versatile numerical tool for developing experimentally validated models that can simulate complex interactions at cellular scale. Using this numerical tool, we presented spring-network models of breast cancer cells that can accurately replicate the experimental data of deformation behavior of the cells flowing in a fluidic domain and passing narrow constrictions comparable to microcapillary. First, using high-speed imaging, we experimentally studied the deformability of breast cancer cell lines with varying metastatic potential (MCF-7 (less invasive), SKBR-3 (medium-high invasive), and MDA-MB-231 (highly invasive)) in terms of their entry time to a constricted microfluidic channel. We observed that MDA-MB-231, that has the highest metastatic potential, is the most deformable cell among the three. Then, by focusing on this cell line, experimental measurements were expanded to two more constricted microchannel dimensions. The experimental deformability data in three constricted microchannel sizes for various cell sizes, enabled accurate identification of the unknown parameters of the spring-network model of the breast cancer cell line (MDA-MB-231). Our results show that the identified parameters depend on the cell size, suggesting the need for a systematic procedure for identifying the size-dependent parameters of spring-network models of cells. As the numerical results show, the presented cell models can simulate the entry process of the cell into constricted channels with very good agreements with the measured experimental data.


Introduction
In addition to several genetic and biochemical factors contributing to the metastatic progression, the biomechanical forces are widely reported to be of considerable importance 1,2 .Among these forces, the hemodynamic forces play a critical role in metastasis progression, for the blood is the primary route of transporting the circulating tumor cells (CTCs) throughout the body 3,4 .Besides, CTC's survival, intravascular arrest, and extravasation are significantly impacted by shear stresses exerted on them originating from their synchronous interactions with the blood plasma, blood cells, and endothelial cells that form the inner layer of blood vessels 1,[5][6][7] .Despite the crucial importance of the impact of such interactions in the hematogenous spread of metastasis, our understanding is far from complete due to the complicated nature of the underlying phenomena.Advanced numerical methods, capable of accounting for the high deformability attribute of CTCs, cell-cell interactions, and interactions with the fluid flow in the delicate microcapillaries have been shown to be promising to decipher CTC's responses to the exerted mechanical forces [8][9][10][11] .Priceless knowledge could be acquired by numerically investigating the impact of various sources of mechanical forces on CTCs.Such knowledge helps answer the perplexing question of when, where, and how a secondary tumor can be formed in a distant organ from migrated CTCs 5,[12][13][14] .Furthermore, it enables the development of more accurate diagnostic methods for localizing secondary tumors and the discovery of more efficient treatments for preventing or delaying this deadly disease 12 .
However, for being trustable, a numerical model should numerically replicate the experimentally observed behavior of CTCs 15 .The deformation behavior of CTCs, that greatly affects every step of the metastatic cascade, has been documented to be one of the key characteristics of malignant CTCs in many research studies 1,[16][17][18] .Thus, the deformation behavior of CTCs can be considered as one of the essential behaviors of CTCs that the numerical model must simulate accurately in comparison with the measured experimental data.Besides, numerical modeling of more complex phenomena such as cell-cell interactions, cell adhesion to vessel walls, and extravasation will be affected by the accurate representation of the deformation behavior of CTCs 9,19,20 .This necessitates developing an experimentally validated numerical model of single CTCs' deformation behavior before moving toward modeling more complex phenomena.This requires acquiring precise measurements of single CTC's deformability in a context similar to their microenvironment.
In vivo models that capture the motion and deformation of CTCs in real time are the most relevant models in terms of replicating the microenvironment 5,21 .Although, challenges with in vivo models such as capturing highresolution images and having precise control of the fluid flow make in vitro models good alternatives for investigating the deformability of cancer cells 12,22 .Since microfluidic devices are the most relevant models in terms of similarity to the blood capillaries among in vitro models, developing advanced microfluidic devices to quantify migratory and deformability capabilities of human cancer cells has attracted many researchers [23][24][25][26][27][28][29] .In this regard, shear flow deformability cytometry 30,31 and extensional flow deformability cytometry 32 are the two contact-free approaches that utilize hydrodynamic flow to change the shape of the cells and assess the cell deformability from cell images captured during cell shape change period 33 .Although these contactless cytometry methods can evaluate cell deformability in a high throughput manner, they are not able to investigate the role of cell-wall interactions on cell deformability.In addition, constricted microfluidic channels whose width are comparable to microcapillary diameters (≈10 μm) 34 became one of the most relevant tools to study the deformability of cancer cells.Using such devices has demonstrated that the entry time, which is the time it takes for the cell to squeeze and fully enter the constriction, would be impacted by the deformability of the cell and is sensitive enough to be different in cancer cells with different levels of malignancy 16,17,28,35 .In addition, since the constricted channel geometry, fluid flow rate, and cells motion are controlled accurately in this type of experiment, it is a good candidate that based on the numerical models of CTC's deformability can be developed and validated.Furthermore, the numerical model is valid when cells are in close contact with the walls of the fluidic domain which could be the case for many practical biomedical applications 21 .More specifically, the cells' entry times calculated from the numerical model can be validated against the experimentally measured ones 8,[36][37][38] .
To date, various numerical methods have been established to quantitatively model the entry process of deformable objects into narrow confinements 37,[39][40][41][42][43][44][45][46] .Most precedent models are continuum models which span from liquid drop models to solid models 47 that considered a cell as a viscous fluid encapsulated with the cell's membrane which obeys related constitutive equations 37,[39][40][41] .However, solid models suffer from accurately modeling large deformation of cells, and liquid drop models are insufficient in simulating cell-wall interactions since they assume a lubrication layer between the cell and the constriction walls 41,[47][48][49] .Although solid continuum models with hyperelastic constitutive laws can improve the accuracy of results for large elastic deformation of cells 50,51 , in general, the applicability of continuum cell models for modeling complex biological phenomena have been limited due to their computational inefficiency in dealing with multiple deformable bodies and integration of biochemical interactions at the cellular scale 52 .
To obviate the mentioned drawbacks, spring network models (SNM) of the cell membrane have been devised, progressed, and applied [53][54][55][56][57][58] .SNM, discretized the cell membrane by a triangular spring network, was devised based on the idea of imitating the physics of the cell membrane by using springs in the numerical model instead of spectrin on the membrane of a red blood cell (RBC) 53,54,59 .High computational efficiency without dropping the model accuracy has been reached with the more advanced SNMs by decreasing the vertices on the cell membrane 55,60,61 .This makes such a method practical for scaling up for simulating circulation, entrapment, and extravasation of CTCs for vessels and organ scale models of metastasis.Moreover, the computational efficiency and accuracy of SNM over continuum models of RBCs suspended in the fluid was demonestrated 49 .More specifically, a comparison between the computational performance of a spring network model and a continuum linear elastic solid model 62 of RBC proved the higher computational efficiency and higher accuracy of the spring network model over the continuum linear solid model.
However, SNMs have a few unknown model parameters, whose quantities affect the deformation behavior of cells directly, need to be identified for each cell type, cell size, and mesh arrangement to replicate the correspondence experiment accurately.The first efforts of finding these parameters were performed on RBC by either manual adjustment or by making a link between SNM and continuum models 57,63,64 .Since these efforts suffered from unrealistic assumptions in identifying the parameters of RBC, the results were not accurate in replicating the experimental data and other rounds of calibration were required 56,64 .In addition, being highly uncertain in finding the unknown parameters of SNM hinders applying this model to other cell types such as CTCs.Recently, a systematic approach called SNM-GA for identifying SNM parameters has been developed to minimize the error between the experimental data and computational calculation using a genetic algorithm (GA) running on supercomputers 8 .SNM-GA is an inverse method that requires appropriate experimental data that should make the identification step feasible and assure producing a validated numerical model able to simulate the deformation behavior of CTCs robustly.Considering the cells' entry process, the data must be reported as cell's entry time and the quantity must be in the order of milliseconds to help SNM-GA identify the parameters feasibly without needing huge computational resources.Besides, measuring the entry time in various devices with various constriction widths increase the number of inputs and outputs that will be used in the parameter identification step and is the key for developing a robust numerical model 8 .To the best of our knowledge the literature lacks such experimental data that reported the entry time of breast cancer cells in order of millisecond using various micro constriction widths.Therefore, capturing such experimental data was one of the goals of this study.
In this study, we first experimentally investigated the deformability of three different breast cancer cell lines (MCF-7, SK-BR-3, and MDA-MB-231) using a constricted microfluidic device whose constriction width is 10 μm.Then, by focusing on the more deformable one (MDA-MB-231), we extended the experiments to two additional constricted microfluidic devices with 8 μm, and 12 μm widths of constriction.Using SNM, the numerical models of a single cancer cell passing through constricted channels have been developed.Afterward, the unknown parameters have been identified based on the measured entry times with the use of GA by minimizing the error between computational calculations and the experimental measurements.We demonstrated that by applying the identified model parameters, both the gradual squeezing of the cell into the constricted channels and the shape change of the cell during the entry process could be accurately replicated in-silico.This study is the first study to present a valid discrete numerical model for the deformation behavior of highly metastatic cancer cells (MDA-MB-231) for a range of cell diameters.The presented models can be used as the foundation of discrete models of the motion and deformation of circulating breast cancer tumor cells upon which numerical models of more complex phenomena such as cell-cell interactions, and CTC cluster deformation behavior in microcapillaries can be built.

Experimental measurement
In this study, to investigate the deformation behavior of breast cancer cells in a microcapillary similar to that of CTCs passing confined spaces, we measured the entry process of the cells entering narrow constricted microfluidic channels.Samples of captured images of the entry process of cancer cells (MDA-MB-231) in all three microfluidic devices (Device #1 to Device #3 from left to right) at five instances are shown in Fig. 1.As shown in Fig. 1b, the entry time starts when the cell starts entering the constriction and ends when the cell fully enters the constriction.In addition, the elongation index, is defined as the ratio of the cell length after completing the entry to the constriction to the original length of the cell before the start of the entry process as shown in Fig. 1b.The images were selected to show similar location and deformation of the cells, but the timing of the squeezing process is different in each device as detailed below.
First, we compared the cell deformation behavior of different breast cancer cell lines MCF-7, SK-BR-3, and MDA-MB-231 with moderate to high metastatic potential 65 (Fig. 2).In Fig. 2a, each point represents one measured cell.As measured data in Fig. 2a show the entry times of all cell lines are increasing exponentially as the cell diameter increases.A comparison between the exponential rate of the best-fit curves of all three cell lines discloses that the entry time of MDA-MB-231 is less sensitive to the increase of cell diameter than the other two cell lines.In addition, entry times of highly invasive breast cancer cells (MDA-MB-231) are shorter than the other two cell lines for a given cell size as the exponential fits in Fig. 2a shows.Taking the cell diameter of 18 μm, the entry times of MDA-MB-231, SK-BR-3, and MCF-7 according to exponential fits in Fig. 2a are 10.6 ms, 16.2 ms, and 27.4 ms, respectively.Figure 2b shows the average entry time for all cell sizes is also shorter for the MDA-MB-231 cell than it is for the other two cell lines.This confirms that the entry time can be a good index to measure the mechanical deformability of cancer cells with various levels of invasiveness 13,14 .Figure 2c illustrates the average entry time for various cell sizes (ranging from 13-21 μm) of the three cell lines.This figure indicates the difference in the deformability of the three cell lines becomes visible for cell sizes greater than 17 μm.The results in Fig. 2 demonstrate higher deformability of MDA-MB-231 cell line compared to the other two cell lines consistent with the studies that reported it to be capable of metastasizing in vivo when it is directly injected to the circulatory system 57 .
Because the deformation behavior of MDA-MB-231 cells entering a narrow channel is closer to that of CTCs passing confined microcapillaries than the other two cell lines, we focused on expanding the deformability measurements of this cell line to also use the data for developing validated SNMs of breast cancer cells.More specifically, we measured the entry time and elongation index of MDA-MB-231 passing through the constricted channel of all three microfluidic devices fabricated in this study.Figure 3 shows the measured entry time for different cell sizes as well as exponential fits for the three devices.The data in Fig. 3 show that for a fixed cell diameter, the entry times and the elongation indexes of the cell passing a wider constricted channel are smaller.For instance, the fitted curves in Fig. 3a show the entry time for an 18 μm MDA-MB-231 that was measured using Device #1, Device#2, Device#3 is 15.1, 10.6, and 6.4 ms, respectively.In addition, the elongation indexes of 18 μm MDA-MB-231, calculated from the fitted curves in Fig. 3b, are 1.27, 1.18, and 1.11 using Device #1, Device #2, and Device #3, respectively.

Parameter identification and numerical results
The deformability of the cells can be captured using the SNM if the correct model parameters are identified.Using the GA, as described in Materials and Methods, we identified the unknown model parameters utilizing the fitted curves to the measured entry times of MDA-MB-231 for 6 different cell sizes between 13 μm to 18 μm in the three microchannel devices (in Fig. 3).In this study, for each cell size the identification procedure was performed using the experimental data of that cell size separately to increase the identification accuracy.For the MDA-MB-231 cancer cell with the diameter of 16 μm, the evolution of the error function (Eq.11) considering all three microchannel devices for the best set of parameters at every generation is shown in Fig. 4a.The numerically calculated entry times in each of the three microchannel devices using this GA 3 µm (passing device #2) and 17.9 µm (passing device #3) at five instances using the three constricted microfluidic devices (Supplementary information Video 1) Entry time starts when cell starts entering the constriction and ends when it fully enters the constriction ðEntry time ¼ t end À t start Þ.Elongation index is the ratio of the cell length after complete entry to the undeformed cell length at the start of the entry process are shown in Fig. 4b.As Fig. 4a shows, for this cell size, the error got close to an optimum value after four generations, improved between generations 5-12, and remained unchanged after the 12th generation.Figure 4b shows that convergence of the entry times follows the same trend as the error function in Fig. 4a. Figure 4c-e shows five instances of the simulated entry process of the 16 μm cell after applying the identified quantities for the unknown parameters of the cell model in three devices.The time instances in these panels are different but were selected to show comparable squeezing state between the three devices.The shape of the cell and its interaction with fluid and microfluidic walls are illustrated from two different views as the cell enters the constricted channel in each microfluidic device.The streamlines in these figures visualize the fluid flow interactions with the squeezing cell and the walls of the devices.Streamlines experience the greatest disturbance in Device #1, where the greatest size mismatch between the cell and constricted channel sizes exists.The disturbance is the most visible in the two last time instances (two rightmost) of Fig. 4c when the cell is near to complete entry to the constricted channel.The same cell size causes a minimal disturbance in the streamlines in Device #3 close to the complete entry state (two leftmost states in Fig. 4e).Supplementary Videos 2 to 4 show the squeezing process of Fig. 4c-e  generation for 16 μm cells, are given in Table 1.Moreover, the results of the identification step including the identified parameter values and the minimized error of every cell size are summarized in Table 2.
The model should be able to replicate the motion and deformation of the population of cells observed experimentally.To verify this, we quantified the numerically calculated entry time and elongation index of the cell after complete entry and compared them with ones obtained experimentally for the population of MDA-MB-231 cells in all three devices represented by the fitted curves in Fig. 5 (reproduced from Fig. 3).
The numerical calculation used the identified parameters reported in Table 2. Figure 5 shows that the numerically obtained entry times and the shape of the cells (characterized by the elongation index at full entry to the constriction) closely replicate the motion and deformation of the population of the cell observed experimentally.The average errors among all devices and cell sizes were 13% for the entry time and 1.4% for the elongation index at full entry to the constriction.
Moreover, in Fig. 6, the numerical and experimental results of the axial cell position and evolution of the cell length during the entry process of MDA-MB-231 in Device #2 were compared.In this figure, data for cell diameters of 13, 15 and 18 µm are presented.Figure 6a shows that increasing the cell diameter makes the cell squeezing into the constriction longer and this behavior was consistently observed in experimental and numerical data.Furthermore, the validated numerical model closely reproduced the experimental observations.Figure 6b shows the gradual and consistent increases in the cell length calculated from the numerical model are very similar to that of the experiment for all three cell sizes.In addition, these observations can also be seen in similarly, Supplementary Videos 5-7 show squeezing of cancer cells into the microchannel, increasingly decelerates the cell motion with increasing the diameter from 13 to 18 µm in the same constricted channel (Device #2). Figure 6c illustrates the axial cell position and cell length at 4 different time steps during cell entry.
We used this validated model to obtain the deformation behavior of the cancer cell in terms of cell length and cell center trajectory using all identified cell model with different sizes (Table 2) that are numerically investigated in detail in Fig. 7.This figure shows that increasing the cell size within one device increasingly decelerates the cell squeezing making the entry time increasingly longer.As an example, in Device #1, the entry time goes up from 4 ms to 10 ms when the cell diameter changes from 13 µm to 18 µm.On the other hand, increasing enlargement of the constricted channel from Device #1 to Device #3 (from 8 µm to 12 µm) increasingly shortens the entry time from the maximum 10 ms for 18 µm cell in Device #1 to ~6 ms for the same cell size in Device #3. Figure 7b shows that the most drastic cell-shape change (5 µm increase in the cell length) occurs with the most cell-channel size mismatch (18-µm-diameter cell in 8-µm channel of Device #1).Consistently, the least cell-shape change (0.5 µm increase in the cell length) occurs with the least cell-channel size mismatch (13-µm-diameter cell in 12-µm channel of Device #3).Within one channel size, increasing the cell diameter increases the cell deformation.For example, in Device #1, cells of diameters 13 µm and 18 µm experience 3.5 µm and 5.3 µm length change, respectively.In addition, using the validated 3D cell models, the effect of changes in the channel height on the entry time of cells with different sizes was investigated (Fig. 8).As Fig. 8 illustrates, for the same fluid flow rate (20 µL/h) flowing through the 10 µm width constriction, increasing the channel height decreases the fluid velocity in the channel, and as a result the cell entry time nonlinearly increases for all three cell sizes.As an example, for the 16 µm cell, increasing the channel height from 22 µm to 30 µm caused the entry time to increase nonlinearly from 2.4 ms to 6.9 ms. Figure 8b shows the side view of the entry process of 18 µm cancer cell at three instances (start, middle and end of cell entry) entering the constricted channel of device #2 with various channel height.
The physical time of these shown instances were different, but all are captured at the comparable instant of the entry process.As this figure shows, the decrease in the fluid velocity magnitude due to increasing the channel height can be noticed specially at the back and front of the cell that illustrates the decrease in cell velocity.
Since in the present model the FSI is a two-way communication between the cell and fluid at each time step 58 , the model is capable of calculating the time-dependent effect of the cell entry process on the fluid flow rate in the microchannel.In fact, capturing the effects of the fluid motion on the cell deformation can help unraveling the impact of hemodynamics on the hematogenous spread of metastasis 6 .Figure 9a shows the changes in the flow rate as a cell squeezes in the constricted channels of the three devices for three cell sizes.In all devices and cell sizes shown in Fig. 9a, the initial flow rate of 20 µL/h starts to drastically drop as cell squeezing starts and progresses but the rate and magnitude of the drop depends on the cell

Conclusion
This study investigated the feasibility of numerically replicating the deformation behavior of single CTCs passing confined spaces in microcapillaries using an advanced in-silico method at the cellular scale.The proposed SNM-based in-silico method was enabled by the previously developed optimization code that encapsulates the experimental data in a validated model of cancer cell deformation.For simplicity, this study used experimental deformation of cancer cell lines in vitro in constricted microfluidic devices instead of studying real CTC deformations in real human microcapillaries.In the first step of this study, the deformability of three different breast cancer cell lines (MCF-7, SK-BR-3, and MDA-MB-231) was measured.The results show the highly metastatic cell line (MDA-MB-231) is the most deformable one among the three confirming this cell line merit for using as CTCs' replica in the in-vitro experiments.Therefore, the optimization algorithm for parameter identification of the cell membrane model was applied to the experimental results of the deformability of highly metastatic breast cancer cells (MDA-Mb-231) passing constricted microfluidic devices with various widths of the constricted channel (ranging from 8 to 12 μm).The parameter identification step, which was performed on various cell sizes (ranging from 13 to 18 μm), helps achieve accurate results for the deformation behavior of discrete cell models in a range of cell diameters.The numerical results show good agreements with the experimental ones in terms of both the entry time and elongation index in various geometrical domains.This means that the results of the numerical models are valid for both the gradual squeezing of the cell into the constrictions and the shape of the cell during the entry process.To the best of our knowledge, this study is the first one to present a valid discrete

Cancer cells culture
Immortalized breast cancer cell vials (MCF-7, SK-BR-3, and MDA-MB-231 cell lines with moderate to high metastatic potential 65 ) were taken out from storage and cultured in Dulbecco's Modified Eagle Medium, DMEM, (4.5 g/L glucose, with L-glutamine & phenol red without sodium pyruvate, WISENT Inc.) supplemented with 10% v/v of fetal bovine serum (WISENT Inc.) and 1% v/v of streptomycin (WISENT Inc.) in 75 cm 2 T-flasks (Thermo Fisher Scientific).Cells were incubated in a standard humidified incubator at 37 °C and 5% CO2.Cells upon reaching a confluency of more than 70% were detached and passaged.On the day of the experiment (half an hour before the experiments), cells were detached with 0.25% trypsin-EDTA (WISENT Inc.) after three or less passages.Cells were quantified for their viability with trypan blue with an automated cell counter (Olympus Life Science).Two to four million cells per mL were grown and used for every experiment.

Design, fabrication, and characterization of microfluidic devices
The design of microfluidic devices contains a single constricted channel in the middle with sizes comparable to human microcapillaries to assure cell deformation at the constrictions' entrance (Fig. 10a, b).The microfluidic devices were fabricated using soft lithography techniques.One layer master mold was fabricated by spin coating SU-8 photoresist 2015 (Kayaku Advanced Materials) at 1400 rpm for 30 s on a silicon wafer.After spin-coating, the wafer was soft baked at 65 °C for 1 min followed by 95 °C for 4 min, and 65 °C for 1 min.Then, the photoresist was exposed to UV light using Karl Suss MA6 Mask Aligner through a previously designed and fabricated chromium glass mask (Nanofab, Alberta University) for 6 s.Next, the post-exposure bake was performed similar to the soft bake procedure.Then, the wafer was developed for 5 min with a SU-8 developer and located on a hot plate at 150 °C for 20 min to stabilize the SU-8 microstructures.The microscopic images of the constricted channels patterned in the master mold are shown in Fig. 10c.In addition, the height of the microchannels was measured to be 28 μm using a surface profilometer (Dektak 8 Stylus Profilometer).After fabricating the silicon master mold, polydimethylsiloxane (PDMS) monomer and curing agent were mixed at a 10:1 volumetric ratio.Then, the mixture was degasified in the desiccator, poured on the silicon master, and thermally cured at 70 °C for 2 h.The cured PDMS was stripped off from the silicon master.Then, the cured PDMS containing the constricted microchannel replicas were cut to the appropriate size, punched in the inlets and outlets to obtain the fluidic access holes, and bounded to a glass slide using oxygen-plasma bonding.The device fabrication was completed by connecting silicon tubing secured with glue to the fluidic access holes.To characterize the critical features of the microfluidic devices, bright field images of the constricted channel of all devices captured by an inverted microscope (Nikon Ti Eclipse) were analyzed.Using a 20× objective, magnified images of the constricted channels, constriction entrance, and exit portion were captured and measured.As Fig. 10a, d show three constricted devices each of which has a 45°t apered entrance at the constricted channel, whose width is comparable to microcapillary diameters ranging from 8 to 12 µm 34 , that fabricated and used in this study.Moreover, comparing to the constricted channels used in the literature 16,17,38 , measuring cell deformability with cell diameter ranging from 13 µm to 26 µm, the fabricated constricted channels in the present study made measuring the single-cell deformability of almost all sizes of the targeted cancer cells feasible.In each device, the width of the channel before the constricted channel is 3 times larger than the width of the constricted channel.This assures both the device fabrication with accurate features and keeping the numerical domain small enough for performing the parameters identification.It is worth mentioning that cell samples were not filtered prior to infusing into the microfluidic devices, but square shape filters with 50 µm distance from each other were devised in the inlet of each device to help reducing cell aggregate at the constricted channel (Fig. 10b).

Numerical method
Lattice Boltzmann method (LBM) The entry process of a single cancer cell into the constricted channels was modeled using Hemocell opensource code (version 2.4) 57 .In this code, the fluid is considered as an incompressible Newtonian fluid whose motion is described in a Eulerian framework and solved by Lattice Boltzmann Method (LBM) implemented in Palabos open-source code (version 2.0) 66 .More specifically, a three-dimensional 19-velocity cube lattice scheme (D3Q19) is utilized in the LBM governing equations as follows: where n i x; t ð Þ, e i , Δt, τ and f i x; t ð Þ are the density distribution function, the velocity vector, time step, relaxation time toward the equilibrium distribution n eq i , and external force, respectively.At each lattice site, the macroscopic fluid density ρ and velocity u can be obtained from the particle density functions as follows: The numerical domains of the three constricted microfluidic devices (shown in Fig. 10a, d) have been created and meshed in a CAD software (Salome 9.7.0)) 67nd the fluid passing through each microchannel was modeled using the above described LBM implementation in Hemocell and Palabos.Figure 10e shows fluid flow simulation at the mid-plane along height direction (z-axis) in three microfluidic devices.Moreover, both the channels' geometry and the geometrical domain in x-y plane that used for numerical simulation of the cell entry process in the constricted devices were shown in Fig. 10e.

Cell Membrane Model
The cancer cell is considered as a membrane with a spherical shape which is discretized by two-dimensional triangles with springs on the triangles' edges.The constitutive equations governing the deformation behavior of the cancer cell include a set of forces (the link force, the bending force, the local area force, and the volume force) acting on the cell membrane as described below 57 .
The link force acts along the edge that connects two adjacent cell membrane vertices is representing the stretch force on the vertices as defined below: where K l , k B , T, L i , L 0 are the link modulus, the Boltzmann constant, temperature, the current length of the edge, and the initial length of the edge, respectively.P ¼ 7:5 nm is the persistence-length of the edge, and τ l ¼ 3 is the relative expansion ratio at which the edge reaches its persistence length 57 .
The bending force is defined in terms of the change in the angles between the normal vectors of two adjacent surface elements as follows: where K b , θ i , and θ 0 are the bending modulus, current and initial angles between the normal vectors of the surface elements, respectively.τ b is the limiting angle and is chosen to be π 6 to prevent unrealistic sharp surface edges 57 .
The local area force applies on each surface element vertices and represents the reaction of the element to change of its area as follows: where K a , A i , and A 0 are the area modulus, the current and the initial area of the triangle, respectively.τ a ¼ 0:3 is the area limiting factor to prohibit surface area changes more than 30% 57 .
The global volume force applies on all vertices of the cell membrane and conserves the volume of the cell.
where K v , V i , and V 0 are the volume modulus the current and the initial volume of the cell membrane, respectively.τ v ¼ 0:01 is the volume limiting factor to resist changes in the cell volume 57 .
It worths noting that the cell model used in this study consists of 642 nodes on which all the mentioned forces have been applied at every time step.Figure 11 illustrates the magnitude of the forces at three instances of cell entry process (start, middle, end) for 18 µm cancer cell model entering the constricted channel of device #2.This figure shows all forces are at their minimum value before the cell deformation starts.For the most or all the vertices the values of the Volume forces, the Link forces, and the Area forces increase as the cell is squeezing to the constriction and reach to their maximum value at the end of the entry process.For the bending forces maximum values were reached during the cell squeezing.Figure 11b depicts the forces act on the cell membrane nodes at the mentioned instances by outputting the cell model during the entry process.
To achieve a realistic numerical model in this study, the internal viscosity of the cell was assumed to be different from the exterior fluid.Therefore, a dimensionless parameter named Viscosity Ratio (VR) was considered as follows:

VR ¼
Interior cell fluid viscosity the fluid viscosity Therefore, K l , K b , K a , K v , and VR are dimensionless parameters that need to be identified accurately to enable the numerical model to replicate the deformation behavior of the cancer cell captured in the correspondence experiment.

Fluid-solid interaction (FSI)
Here, fluid-cell interactions were modeled with Immerse Boundary Method (IBM) which acts as a bridge between Eulerian grids of the fluid and Lagrangian grids of the cell membrane 68 .More specifically, the exerted forces on the cell membrane nodes, determined by the cell's constitutive equations ðFÞ, were spread on the fluid grids as follows: where δ is the Dirac delta function, x is the coordinate of the Eulerian grids, and Xðq; tÞ is the position of a cell node with Lagrangian coordinate q at time t.
The velocity of the cell membrane nodes U X q; t ð Þ ð Þwas obtained from the integral below and applied for updating the positions of the nodes.
where u x; t ð Þ is the velocity of the fluid with Eulerian grid x at time t.

Cell-wall interaction
Furthermore, repulsive forces were defined between the nodes of the cell and microchannel wall to avoid cell penetration into the microchannel wall to model the behavior of the cell near the walls as: where κ rep is the repulsion constant, d is the distance between the nodes of cell and wall, d cut is the threshold of repulsive force activation, and m is the unit vector pointing from the wall node to the cell node 57 .The repulsive-force parameters were constant in all simulations in this study (Table 3).

Genetic algorithm (GA)
The mentioned parameters which describe the deformation behavior of the cancer cell were identified using the previously developed genetic algorithm 8 .This algorithm benefits from creating the first population of size 60 using randomly generated multi-digits binary numbers.Every 128 digits of the binary numbers represents one of the parameters and each row of the first population represents a set of parameters K l , K b , K a , K v , and VR.Afterwards, every two rows of the first population were selected as parents and crossover has been performed for the crossover probability higher than a user-defined probability (0.8) for generating children.Then, mutation step was conducted for the mutation probability higher than a user-defined probability (0.6) on every digit of the binary numbers that made by crossover.The made children and the mutated ones were added to the previously generated first population.Then, the binary numbers were converted to decimal numbers according to the upper and lower bounds of each parameter as provided in Table 4.Then, using the decimal numbers, numerical simulations of the cancer cell entry process were performed in all three devices benefiting from parallel jobs on supercomputers.It is worth noting that every generation consisted of 120 sets of parameters, therefore, 360 simulations were performed simultaneously using 2 cores for each and the entry time was stored when the cells fully entered the constricted channels.For those simulations, that the entry process took much longer than the experimental data, the simulations were stopped when the entry time reached a user-defined value (twice of the experimental data).Finally, the outcomes of the numerical simulations and the experimental data were compared based on the below error function: where E n , ET s , and ET e are the error in the nth device, numerically calculated entry time, and experimentally measured entry time, respectively.Here, n t is equal to 3 for using the entry time at three different constricted devices.
At the end, the results for the parameters of the cell constitutive equations were sorted and the best 20 ones added to the initial population of the next generation.The algorithm stops if the error remains unchanged after 20 successive generations.

Experimental setup
Before running the experiments, the devices were degassed for up to one day with Pluronic solution to avoid cell adhesion to the channels and washed with a constant flow of PBS for 20 min.The microfluidic device was placed on the stage of a Nikon Ti Eclipse inverted microscope in the experimental setup shown in Fig. 1a.The cell sample flowing in the media consisting of RPMI-1640 solution with 20% fetal bovine serum (FBS) (the media density ρ ¼ 1020 ± 5 kg m 3 , and the media dynamic viscosity μ ¼ 1:089 ± 0:044 mPa s) 69 with cell viability of % 90 or higher was infused into the constricted channel using a syringe pump (Chemyx Inc., USA) at a constant flow rate of 20 μL/h.Cell samples with the concentration of 2 × 10 6 cells per mL of media were prepared and used for all experiments.In addition, since the main focus of the present study is to validate the numerical model of single cancer cell deformation behavior, only the data of single cancer cells that pass the constricted channels one at a time have been gathered.More specifically, the captured data of cell clusters and more than one cells in the constricted channels at the same time were set aside from further analysis.At this flow rate, the entry time measured for the average cell size is less than 16 ms.The flow of  cancer cells into the constricted channels has visualized in bright field mode of the inverted microscope using a ×20 magnification objective and recorded using a high-speed camera (FASTCAM S1 model, Photon USA, Inc.) at a high frame rate of 5000 fps with the spatial resolution of 512 × 512 pixels.All captured videos and images were analyzed manually using Photron Fastcam Viewer 4 (PFV4) software to measure cell size, entry time, and elongation index which is the ratio of cell length after entering the constriction to the original cell length.

Fig. 1
Fig. 1 Experimental set up for cancer cell deformability measurement.a Experimental set up for high-speed measurements of the single breast cancer cell deformability.b Captured images from entry process of the single breast cancer cell with the diameter of 18.2 µm (passing device #1), 22.3 µm (passing device #2) and 17.9 µm (passing device #3) at five instances using the three constricted microfluidic devices (Supplementary information Video 1) Entry time starts when cell starts entering the constriction and ends when it fully enters the constriction ðEntry time ¼ t end À t start Þ.Elongation index is the ratio of the cell length after complete entry to the undeformed cell length at the start of the entry ðElongation index ¼ L1 L0 Þ

Fig. 2
Fig. 2 Deformability measurement of breast cancer cell lines with various level of invasiveness.Comparison of the deformability of breast cancer cell lines (MCF-7, SK-BR-3, and MDA-MB-231) using device #2 by illustrating the measured data with a scatter graph showing all points of entry time vs cell diameter for every single cell, b box plot of entry time of all captured data in the three cell lines, and c the heat map for average entry time in various cell diameter

Fig. 3 31 Fig. 4
Fig. 3 Deformability measurement of MDA-MB-231 using the three microfluidic devices.Experimental measurement of the deformation behavior of the highly metastatic breast cancer cell line (MDA-MB-231) using the three microfluidic devices by measuring a the entry time of single cancer cell passing the constricted channels, and b the elongation index and microchannel sizes.For example, for the cell diameter of 18 μm in Device #1 the flow rate decreases to 0.7 μL/h while the lowest flow rate in this device for 16 μm cell and 14 μm cell are 3.3 μL/h and 7.4, respectively.Besides, the constricted channel width affects the fluid flow rate decrease during the entry process as for 18 μm cell in Device #2 and Device #3 the flow rate reaches as low as 3.8 μL/h, and 5.8 μL/h, respectively.Figure9bshows the effect of the flow rate on the entry time of the cell.For the flow rates more than 40 µl/h, the entry times of three different cell sizes (14, 16, 18 µm) are almost the same in the constricted channel of all three devices.For instance, the numerically calculated entry times for cell diameters of 14, 16, and 18 μm that pass the constriction in Device #1 with 40 μL/h flow rate are 1.4 ms, 1.8 ms, and 2.3 ms, respectively.As the flow rate decreased from 40 μL/h to lower values, the effects of the cell and microchannel sizes become increasingly more remarkable.For example, in Device #1 with the flow rate of 20 μL/h the entry times of 4.2, 7.3, and 9.6 ms

Fig. 5
Fig. 5 Comparison between numerical results with the curve fitted on the population of the experimental data.The comparison between experimental results and the numerical ones for a the entry time, and b elongation index of MDA-MB-231 cell line taken after applying the identified parameters of the cell model

9 Fig. 6 Fig. 7
Fig. 6 Comparison between numerical results and experimental ones for cell deformation behavior of three distinct data.The comparison between experimental results and the numerical ones for a axial position, and b cell length of three MDA-MB-231 cells with different size during their entry into constricted channel of device #2.(Supplementary Video 5 to Video 7) c Descriptive images of cell length change and cell axial position during the entry process

Fig. 8
Fig. 8 Numerical investigation of channel height effect on cells' entry times.a The numerical results of the cancer cells entry times entering the 10 µm constricted channel with various heights at 20 µL/h fluid flow rate using the developed cell models.b The side view images of the numerically simulated entry process of the cell entering the 10 µm constricted channels with the channel height 22, 26, 30 µm at start, middle, and end of the entry process

Fig. 9 Fig. 10
Fig. 9 Numerical investigation of fluid flow rate change during cell entry process and the effect of fluid flow rate on cells' entry time.a Change of the flow rate while the cell entering the channel that was calculated numerically indicates flow rate change dependency to the cell size.b The effect of fluid flow rate on the cell entry time of various cell sizes in the three different constricted channels was investigated numerically

Fig. 11
Fig. 11 Magnitude of forces acting on the cell membrane during cell entry.a Magnitude of cell model forces acting on every node of the cell extracted at three different instances (start, middle, end) during cell entry process for the 18 µm cell entering the constricted channel of device #2.b Descriptive images of the cell model forces contributing to entrance of the cell into the constriction at the same instances as part (a) Figure 1b-d shows MDA-MB-231 cell squeezing and entering the microchannel sizes of 8, 10 and 12 μm, respectively, at five instances until completion of cancer cell entry.

Table 1
The quantities of the 16 μm cell model parameters for the best ones at each generation

Table 2
The quantities of the identified parameters for various cell sizes ranging from(13-18 μm)

Table 3
Parameters used in this study

Table 4
Upper bound and lower bound for cell model parameters